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This paper presents a detailed, numerical study on the performance of the 
^| standard phasing algorithms with random phase illumination (RPI). Phasing 

^ ' with high resolution RPI and the oversampling ratio a = 4 determines a unique 

phasing solution up to a global phase factor. Under this condition, the standard 
phasing algorithms converge rapidly to the true solution without stagnation. 
Excellent approximation is achieved after a small number of iterations, not just 
^ with high resolution but also low resolution RPI in the presence of additive 

as well multiplicative noises. It is shown that RPI with a = 2 is sufficient for 
^ \ phasing complex-valued images under a sector condition and a = 1 for phasing 

nonnegative images. The Error Reduction algorithm with RPI is proved to 
converge to the true solution under proper conditions. 
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1. Introduction 

Fourier phase retrieval is the problem of reconstructing an unknown image from its Fourier 
magnitude data. Phase retrieval is fundamental in many applications such as X-ray crystal- 
lography [2], astronomy [3], coherent light microscopy [3], quantum state tomography and 
remote sensing. 
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Due to the absence of the phase information, phase retrieval does not have a unique 
solution. Phase retrieval literature has long settled with the notion of uniqueness modulo 
the trivial ambiguities of spatial shift, conjugate inversion and global phase [5] [6] and focused 
on circumventing the stagnation problem associated with the standard phasing algorithms. 
The numerical stagnation problem is often attributed to the nonconvex constraint imposed 
by the Fourier magnitude data [7] [S] [H] [ID] ■ 

In this paper, we explore a phasing method based on random phase modulator which 
randomly modifies the phases of the original image by a mask. As proved in [1] phasing with 
random (phase or amplitude) illumination often leads to a unique solution up to a global 
phase factor (here dubbed absolute uniqueness) . In what follows we show that phasing with 
random phase illumination (RPI) also leads to superior numerical performances, including 
rapid convergence, much reduced data and noise stability of the standard algorithms. We 
show that under proper conditions the Error-Reduction (ER) algorithm with RPI converges 
to the true solution (Theorem H]). 

Consider the discrete version of the phase retrieval problem: Let n = (ni, . . . ,71^) G Z d 
and z = (zi, . . . , Zd) G C d . Define the multi-index notation z n = z^z^ 2 ■ ■ ■ z n d d ■ Let C(N) 
denote the set of finite complex-valued functions on Z d vanishing outside 

M = {0 < n < N}, N = (N u N 2 ,..., N d ). 

d 

Here m < n if rrij < rij,\/j. Denote \J\f\ = J~J-Wj- 

3=1 

The ^-transform of a d dimensional finite array /(n) G C(N) is given by 

F(z) = ^/(n)z- n . 

n 

The Fourier transform can be obtained from the ^-transform as 

F(e i2 ™) = /(n)e~ 2 ™ (1) 

n 

for uj = (ui, oj 2 , ■ ■ ■ , uJd), < Uj < 1. 
From the calculation 

N 

\F(e*™)\ 2 = E E /(m + n)J(^)e- i2 ™ w 

n=-N m+neAT 

we see that the Fourier magnitude measurement is equivalent to the standard discrete Fourier 
measurement of the correlation function 

C/(n)= E/( m + n )7R (2) 

meJV 
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if sampled at the lattice 

f \ i 1 2 2^ -1 

I ' 2A^ + 1 ' 27V,- + 1 2iVj + 1 J V ; 

which is 2 d times of the grid of the original image. The standard phasing problem is to 
recover the array /(n) from its Fourier intensity measurement Y{ui) = \F(e l27TU} )\ for u) £ £ 
or smaller sampling sets. 

Clearly the correlation function Cf and the Fourier magnitude data are invariant under 
spatial translation 

/(•) -> /(• + t) for some t £ Z d , 

conjugate inversion 

/(.)->7(N— ) 

and constant global phase change 

These trivial associates all share the same global geometric information as the original object. 
The classical results of uniqueness given in [5] [6] p2] say that for almost all objects in 
dimension two or higher the trivial associates are the only ambiguities there are with phase 
retrieval. 

On the other hand, by dimension counting Miao et al. [TT] have argued that overall 2 
times oversampling, independent of the dimension d, uniquely determines a unique phasing 
solution up to spatial shift, conjugate inversion and global phase factor. To measure the 
degree of oversampling we use the oversampling ratio (OR) 

Fourier magnitude data number 
unknown-valued image pixel number 

introduced in [TT]. As we demonstrate below, Miao et al.'s conjecture can be realized by 
using RPI, but not uniform illumination. 

As shown in [1] random illumination (RI) can help remove the phasing ambiguities of 
spatial shift and conjugate inversion. An illumination amounts to replacing the original 
image /(n) by 

g(n) = A(n)/(n), 

where A(n) is a known array representing the incident wave. In the case of uniform illumi- 
nation, A(n) = 1. In the case of random phase illumination (RPI) |13j . 

A(n) = e^ (n) (4) 
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where 0(n) are random variables on [0, 27r], and in the case of random amplitude illumi- 
nation [Hl[l5], -M n ) * s an arra y °f rea l random variables. RI can be facilitated by random 
phase/amplitude modulators or random masks. 

The paper is organized as follows. We review the absolute uniqueness of phasing with 
RPI in Section [2] and standard phasing algorithms in Section [3] where convergence of the 
Error Reduction (ER) iteration to the true solution is presented (Theorem H]). We present 
the numerical phasing results in Section |H We conclude in Section 

For the rest of the paper we use the following notation: the vector space C{M) is endowed 
with the inner product < f,g >= Y2 n f( n )9( n )- ^ OT a complex number z, ^ft(z) and Sf(z) 
denote the real and imaginary part of z. Lz G [0, 27r) denotes the phase (angle) of z. When 
z — 0, Lz is taken to be unless specified otherwise, [a] = a(mod(27r)). 

2. Uniqueness 

In the following we recall several uniqueness results from [T] relevant to phasing with RPI. 

First we define the rank of an array. The support of the array consists of the set of nonzero 
pixels. The rank of the array is the dimension of its support's convex hull in M d . 

Theorem 1. Let A(n) be independent, continuous random variables on S . Let /(n) G C{M) 
be a real-valued array of rank > 2. Then, with probability one, f is determined absolutely 
uniquely up to ± sign by the Fourier magnitude measurement on C. 

A more general, practical constraint is to restrict the image values within a certain sector 
of the complex plane. For instance, when the incident X-rays are low energy photons(soft 
X-rays), the electron density is complex. The real part represents the effective number of 
electrons that diffract the X-rays in phase and is usually positive but becomes negative only 
when the energy of the incident X-rays is near an absorption edge. The imaginary part 
represents the absorption of the X-rays by the specimen and thus is always positive. 

Theorem 2. Let A(n) be independent, continuous random variables on S 1 . Let f be a 
complex-valued array of rank > 2 such that Z/(n) G [a,/3],Vn. Let S denote the sparsity of 
the image and let U_5'/2JJ be the greatest integer less than or equal to S/2. 

Suppose that the phases 0(n) of RPI are independent, uniform random variables on [0, 27r]. 
Then with probability no less than 1 — |A/"|(/3 — a)^ s ^ 2 ^(27r)~^- s ^ 2 ^ , the object f is uniquely 
determined, up to a global phase, by the Fourier magnitude measurement on C. 

The global phase is uniquely determined if the angular sector [a, 0\ is tight in the sense 
that no proper subset of [a, b] contains all the phases of the object. 
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For general complex- valued images without any constraint, we use two independent RPIs 
to collect data. 

Theorem 3. Let Ai(n) and A2(n) be two independent arrays of continuous random variables 
on S 1 . Let /(n) G C(Af) be any complex-valued array of rank > 2. Then almost surely /(n) is 
uniquely determined, up to a constant phase factor, by the Fourier magnitude measurement 
on £ with two illuminations Ai and A2. If the second illumination A2(n) is deterministic 
while Ai(n) is random as above, then the same conclusion holds. 

3. Phasing Algorithms 

To find the true object satisfying both the object-domain constraint, which is usually convex, 
and the frequency-domain constraint, which is non-convex, most phasing algorithms are 
based on the idea of alternating projections from the convexity literature [9]. 

3. A. Projections 

Defintion 1. Let V be a subset of C(N), the orthogonal projection of f G C{M) on T> is 
argmin geV \\g - f\\ 

If the minimizer is not unique, one of them is arbitrarily selected. When T> is a closed 
convex subset oiC(J\f), the minimizer is unique. 

Proposition 1. Let T> denote any closed convex subset of C(N) and let f be any element 
in C(N). Then there exits a unique h G X> such that 

inf gev \\9 ~ f\\ = \\h ~ f\l 



Let r be the set of functions satisfying the object-domain constraint, such as a known 
support or positivity, and fl be the set of functions satisfying the frequency- domain constraint 
imposed by the known Fourier magnitude data. A solution of phase retrieval is a function 
belonging to T D Q. Let V and Vf be the orthogonal projection on T and Q respectively. 

Let A be the diagonal matrix with diagonal elements A(n), and set g = Af. Let $ be the 
discrete Fourier transform and set Y = |3>/|. 

Given the Fourier intensity data Y, we define the intensity fitting operator T as 

G\u)=T{G}^) = { if !^)l>° . (5) 

if |G(w)| = v ' 
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When G{(jj) = 0, ZG(u;) is not uniquely defined and ZG(ui) is set in In this case, 
Indeed /CG(u)) can be arbitrarily chosen at the zero set of G, and we define 



Y(u>)e UG W if|G(w)|>0 
Y(u)e i0 W if |G(u>)| = 



where 

T e {G}{u) = 

The object domain projection V Q can take a varied form depending on the problem. 

• When T is the set of images with a given phase a, 

V {h}(n) =V a {h(n)} = max{3(/i(n)) sin a + U(h{n)) cos a, 0}e ia . 

• When T is the set of images with phases in [a, f3] for < a < (3 < 2tt, 



if (3 -a < 7T, V {h}(n) 



' h{n) if a ~< lh{n) -< (3 

Vp{h(n)} if f3 -< lh(n) -< [13 + tt/2] 

M(V a {h{n)}) if [a - tt/2] -< Z/i(n) -< f3 

else 



/i(n) if a -< Z/i(n) -< /3 

if /3 - a > 7T, P Q {/i}(n) = <j Vp{h(n)} if £ -< Z/i(n) -< [(a + /3)/2 + tt] 

P a {/i(n)} if [(a + /3)/2 + tt] -< Z/i(n) -< a 

where a -< 6 -<b means 9 is between a and 6 such that 

a<6 <b iia<b 
a<6 < 2tt or 0<6 <b if a > b 

When T is the set of real valued images, 

V {h}(n) = R(h(n)). 

When T is the set of nonnegative real-valued images, 

V {h}(n) = max{K(/i(n)),0}. 
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When T is the set of complex valued images with nonnegative real and imaginary parts, 

^(V {h}(n)) = max(5ft(/i(n)),0) 

%(V {h}(n)) = max(3(/i(n)),0). 
When T is the set of images with support S, 



V {h}(n) 



h(n) if n G S 
else 



Two error metrics e Q and Sf defined by 

e (h) = \\V {h} - h\\, 

e f (h) = \\V f {h}-h\\ 
play an important role of our studies. When $A is unitary, as in the case of RPI, 

£f (h) = \\V f {h} - h\\ = \\T®Ah - $Ah\\ = \\Y — \$Ah\ \\. 

3.B. Over sampling 

The oversampling method has proven to be an effective, flexible way of implementing var- 
ious phasing algorithms by converting Fourier magnitude data more finely sampled than 
demanded by the original image grid into zero padding which then acts like a support con- 
straint of the original image [SI UHl EOl El] . In this set-up, the oversampling ratio is given 
by 

image pixel number + zero-padding pixel number 

a = . 

image pixel number 

3. C. Error reduction ( ER ) 

ER algorithm [17] is based on the Gerchberg-Saxton algorithm [18] and is the most basic 



phasing algorithm. ER is the plain version of the alternated projection method: 

f k+1 = V V f f k (8) 

which can be conveniently represented by the following diagram 

ER enjoys the error-decreasing property following the same argument in [TT] . 
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fk+l/ fk 



Vo 



fk 



•9k 



■9k 



T 



Fig. 1. 



Proposition 2. Let V be a closed convex subset of C(N). Let $ and A be unitary matrices. 
Then the array {fk} produced in (jSJ) satisfies 



£/(/fe+i) < e/(/*). 
The equality holds if and only if fk+i = fk ■ 
Proof. 



(9) 



> 



> 



II/* — /fell 
||/fe+i - /fell 
||Cfc+i — GfcH 
||Gfc + i — Gfc+ill 
= ||/fe+i — /fc+ill 

= £/(/fc+i). 

The equality holds only if ||/ fe - = ||/ fc+ i -/£||, where = ^o{/fe}- Since T is a closed 
convex subset, = /j, according to Proposition [TJ □ 

Remark 1. Proposition® holds for the fk+i = V Q Vjfk with arbitrary 9(uj). 

Proposition [2] shows that the error £/(/fe) decreases strictly until it reaches a fixed point 
of V Vf, implying that the ER iteration converges to a fixed point. 

Proposition 3. Let fk+i = V Vffk- Let T be a closed convex subset of C{M) and $,A be 
unitary matrices. Then every convergent subsequence of {fk} converges to some h such that 

1. if $A/i(u;) 7^ 0, Vo; E C, h is a fixed point ofV Vf. 
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2. if $A/i(u;) = for some u) G C, h is a fixed point of V Q Vj for some 9. 

The proof of Proposition [3] is given in the Appendix. The question is, Is a fixed point of 
ER necessarily a phasing solution? With the uniform illumination, however, this is generally 
not true [19]. When a fixed point fails to be a phasing solution, it is called a trap and can 
plague the reconstruction procedure (cf. Figure 0(a), H]( a) and|U(c)). 

Below, we answer this question in the affirmative under certain assumptions for the case 
of RPI. The difficulty is ER may converge to a fixed point of V Vj which fails to satisfy the 
Fourier magnitude data. In other words, the limiting point h may not be a fixed point of V %. 

In the following main theoretical result of the paper, we prove that if Vjh satisfies the 
zero-padding condition, then it must be the phasing solution. 

Theorem 4. Let /(n) G C{M) be an array with /(0) 7^ and of rank > 2. Let A(n) be i.i.d. 
continuous random variables on S 1 . Let the Fourier magnitude be sampled on C Let h be a 
fixed point ofV Vj such that Vjh satisfies the zero-padding condition. 

(a) If f is real-valued, h = ±/ with probability one, 

(b) If f satisfies the sector condition of Theorem^ then h = e lu f , for some v, and satisfies 

the same sector constraint with probability at least 1 — \J\f\(/3 — a)^ 5 / 2 "(27r) - ^ 5 / 2 ". 



The hybrid input-output (HIO) algorithm is a widely used, better-performing phasing 
method than ER's [17]. HIO differs from ER in how to update the image in the object 
domain in order to avoid the trapping and stagnation. 

Below we present a modified version of Fienup's HIO which performs better than the 
original version. We refer to Figure [1] for the notation. In HIO, the last step V of ER 
iteration is replaced by the following. 

• When T is the set of real-valued images, 



3.D. HIO 



HC/WiCn)) = K(/i(n)) 



(10) 
(11) 



If, in addition, the nonnegativity constraint is assumed, then 




(12) 



(13) 
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• When T is the set of complex- valued images with nonnegative real and imaginary parts, 

Mf (n)) = I m ' k{n)) if»(/£(n))>0 
Ufe+U >} \ n(f k (n))-P-®(f> k (n)) if < " 1 } 

^ +1 (n))4^ (n)) (15 ) 

U 1 ^ 1 »(/ fc (n))-/3-^(n)) if 9f(/i(n)) < 1 J 

Algorithms with two illuminations 

Let Ai(n) and A2(n) be two arrays representing two illuminating fields. Two sets of Fourier 
magnitude data Y\ = \$A\f\ and Y\ = \$A 2 f\ are collected, each with an OR a. Let T\ and 
T 2 be the intensity fitting operators corresponding to Y\ and Y 2 , respectively, as in Thus 
the projections onto the set of images satisfying the Fourier magnitude data Y\ and Y 2 are, 
respectively, 

V 1 = Ai ^-Vi^Ai 

and 

The corresponding ER algorithm with two sets of Fourier magnitude data Y\ and Y 2 is 
given by 

fk+i = V V 2 VJ k . (16) 
The corresponding HIO is obtained by replacing V in ( Till by ( frol - (fT5l) . 
4. Numerical Simulations 

In this section, we perform numerical phasing from the Fourier intensity measurement with 
UI or RPI. 

Our test images are the 256 x 256 Cameraman and the 138 x 184 Phantom. We surround 
both images by dark (i.e. zero-valued) border to create images of loose support. Images of 
loose support are typically more challenging to reconstruct. For Cameraman the border is 
13 pixel wide in each dimension and the resulting image has 269 x 269 pixels in total. For 
Phantom the dark margin is such that the resulting image has 200 x 200 pixels. 

For the oversampling ratio a, we zero pad the images to generate a 269 y/a x 269\/a 
Cameraman and 200y/a x 200^/cr Phantom. We synthesize the Fourier magnitude data by 
applying the FFT to the array. 
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(a) 



(b) 



Fig. 2. Test images of loose support: (a) 269 x 269 Cameraman (b) 200 x 200 
Phantom where the dark borders represent loose support. 



4- A. Error, Residual and Noise 

Let / be the recovered image. The relative error is defined as 

II/ - /II /II /II ^ absolute uniqueness holds 

if uniqueness holds only up to a global phase 



e(f) 



min \\f-e tv f 



and the relative residual is defined as 

r(f) - 



Y-\$AV {f}\ 
\\Y\\ 



where V Q is introduced if / may not strictly satisfy the object domain constraint as in the 
case of HIO. 

We consider three types of noise: Gaussian, Poisson and illumination noise, the last of 
which is defined as follows. Suppose the illumination field is noisy A(n) = exp(i0(n)) 
with 0(n) = 0(n) + t(5, n) where t(S, n) are independent, uniform random variables in 
[-thJ/IOO, Ti-5/100], 5 > 0. 

We also test phasing with low resolution illumination which does not consist independently 
distributed pixel values but independently distributed blocks of deterministic (indeed, uni- 
form) values. In our experiments, illumination of independent 40 x 40 blocks works well for 
real-valued nonnegative images and, for complex images, illumination of independent 4x4 
blocks works well. 
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4-B. Convergence Test 

The reconstruction of the real- valued nonnegative images Cameraman and Phantom with one 
UI or RPI is shown in Figures[3]andH]respectively. For Figures [3] and HI we terminate the pure 
ERwhen ||/ fc+ i-/ fc ||/||/ fc || < 0.01%. For HIO+ER, HIO is stopped when ||/ fc +i-MI/IIMI < 
1% with a maximal 1000 iterations and ER is terminated when ||/fc+i — /fc||/||/fc|| < 0.01%. 
In Figure E] and HI \\f — V Vff\\/\\f\\ is as small as 0.01%, implying that / is near a fixed 
point of V Vf. 

As commented before the pure ER iteration always converges to a fixed point of V Vf. But 
with one uniform illumination and a = 4, the fixed point of V Vf is not a phasing solution 
as the relative residual stagnates at 5.49% in Figure E(b) and at 14.71% in Figure BJb). HIO 
followed by ER improves the recovery over pure ER but the recovered Cameraman in Figure 
[3](c) displays the well known artifact of stripe pattern and the recovered Phantom in Figure 
S(c) is severely blurred and distored. 

With one low resolution RPI (block size: 40 x 40) and a = 2, the recovered images in 
Figure Ete), E(g), Hte) andH(g) are excellent approximation to the true images, even though 
absolute uniqueness is not guaranteed for low resolution RPI. HIO+ER is superior to pure 
ER in significant speed-up in convergence (Figure E(f) versus E(h), Figure HJ^f) versus E(h)). 

With one high resolution RPI, high quality reconstruction can still be achieved with the 
oversampling ratio equal to 1, cf. Figures E(i), E(k), E(i) and@Jk). Notice the rapid conver- 
gence of HIO+ER in Figures [3(1) and Hp). 

4-C. Oversampling Ratio Test 

To systematically test the oversampling ratio required for phasing with RPI, we introduce 
5% different types of noise (Gaussian, Poisson, Illumination), use low (block size: 40 x 40) 
as well as high resolution RPI and let a vary. We use an adaptive version of HIO+ER: HIO 
and ER are terminated if the residual increases in 5 consecutive iterations. The relative error 
of reconstruction for the nonnegative image Phantom is averaged over 5 trials and shown in 
Figure [5(a). Clearly the relative error steadily decreases as the oversampling ratio increases. 
Without noise, low resolution RPI can achieve near zero error with a = 1.1. With 5% noise, 
the relative error stabilizes after a = 2 to a level comparable to the noise. 

Next we consider the complex-valued Phantom with phases randomly distributed in the 
sector [0, 7r/2]. Figure [5(b) shows the average relative error e(/) with one high resolution 
or low resolution (block size: 4x4) RPI and three kinds of noise. Again the relative error 
stabilizes after a = 2 to a level comparable to the noise. Note that for 1.8 < a < 2, there 
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are more free variables in the complex-valued image than in the Fourier intensity data and 
yet the reconstructions are still of good quality 

Finally, we consider the complex-valued Phantom with phases random distributed in 
[0, 27r]. Figure [5](c) shows the average relative error e(f) with one high resolution or low 
resolution (block size: 4x4) RPI plus one UI. Excellent recovery is achieved for a > 1.8. 

4-D. Stability Test 

For images with positivity constraint and with one RPI, we terminate HIO when the relative 
residual increases for 5 consecutive steps and apply 10 steps of ER afterward. The maximal 
HIO iteration is set to be 100. For complex-valued images with two illuminations, we apply 
200 steps of HIO and 300 steps of ER. 

Figure E] shows the recovery for the nonnegative- valued images with one high resolution 
RPI and 5% Gaussian ((a)-(d)), Poisson ((e)-(h)) and illuminator noise ((i)-(l)). Multiplica- 
tive noise such as Poisson and illumination noises are generally more debilitating than the 
additive Gaussian noise. 

With a low resolution RPI (block size: 40 x 40), the quality of reconstruction suffers slightly 
as shown in Figure [7] for nonnegative- valued images. The deterioration is most visible in the 
case of Poisson noise with the blocky pattern in Figure 0(e) and[7](g). 

Figure [8] shows the average relative error e(f) versus noise for (a) nonnegative-valued 
Phantom and a = 2, (b) Phantom with phases randomly distributed in [0,7r/2] and a = 4 
and (c) Phantom with phases randomly distributed in [0, 27r] and a = 3. One high or low 
(40 x 40) resolution RPI is used in (a) while one high or low (4 x 4) RPI and one UI are 
used in (b) and (c). The adaptive HIO +50 ER is used for (a) and (b) while 200 HIO + 300 
ER is used for (c). 

Relative error increases almost linearly with respect to the relative noise level with the 
noise amplification constant at worst 2. Clearly the illumination noise is most debilitating, 
followed by the Poisson noise. Nevertheless, the noise stability is achieved with even the low 
resolution RPI for all three types of noise. 

5. Conclusion 

We have given a proof of convergence of ER (Theorem H]) and demonstrated that the stag- 
nation problem of standard phasing algorithms such as ER and HIO can be alleviated if the 
ambiguities associated with spatial translation and conjugate inversion are removed by RPI. 
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In addition, phasing with RPI has the following advantages: (i) It is stable with respect to 
additive as well as multiplicative noises with a moderate noise amplification constant; (ii) 
It reduces the oversampling ratio by more than a factor of 2; (iii) It reduces the number of 
iterations by more than an order of magnitude. We have also shown that phasing with RPI 
performs well with low resolution illumination and can tolerate a high level of illumination 
error, adding assurance that the random illumination needs not be calibrated exactly. 

The lower bound a > 2 for phasing of [TT] was never actually achieved but we have achieved 
the lower limit in phasing with RPI for complex- valued images under a sector constraint. For 
nonnegative- valued images, phasing with one high resolution RPI reduces the oversampling 
ratio to unity, the minimum level by the dimensional count. 

Appendices 

A. Proof of Proposition [3] 

Proof. By Proposition |2} lim^oo £f(fk) = Vi f° r some rj > 0. Since fk+i = V {fk}, we have 
||/fe+i|| — WfLW = ll^fcll = ll^ll and that {f k } is a bounded sequence. Every bounded sequence 
in C(Af) has a convergent subsequence, so {f k } has at least one convergent subsequence. 
Without loss of generality, we assume lim^oo f k = f*. Next, we prove that /* must be a 
fixed point of V Vf or V Vj for some 6. 

Since $ and A are unitary matrices, lim^oo $A/fc = $A/* in || • ||, and thus 

lim ®Af k (u) = $A/*(w), Vw. 

• If $A/*(tt>) vanishes nowhere in C, then 

lim Z$A/ fc = Z$A/* 

implying 

lim G' k = lim T$A/ fc = T$A/*. 

k— >oo k— >oo 

Therefore, 

lim f k+1 = V V f f* 

which along with the convergence of £f(fk) and f k implies that £/(/*) = £f(PoPff*)- 
By Proposition |2l we have 

r = v v f r. 
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• If $A/*(tt?) = at some u> G £, T$A/fc(u;) may not converge. However, since T$A/^ is 
bounded in view of ||T$A/fe|| = \\Y\\, there exists a subsequence {fkA and some 9(f*) 
such that limfc^oo 7"$ A/*, (w) = l~ e &Af*(uj) where T e is defined in (jJJ). Therefore 

lim VoA-^^TQAfa. = V Ar l <S>~ l T e <$>kf\ (17) 

namely 

lim = lim V Q V s h 3 = V V]f 

j—toc j—^oo 

which along with the convergence of ^f(fk) an d fk implies that £/(/*) = Ef(V Vjf*). 
By Proposition [21 it follows that 

□ 

B. Proof of Theorem [4] 

Define 

/ m+ (.)=/(m+0, /m-(-) = /(m-0- 

Let 

F(z) = ^/(n)z- n 

n 

be the z-transform of /. According to the fundamental theorem of algebra, F(z) can be 
written uniquely as 

p 

F(z) = az~ no l[F k (z), 

k=l 

where no is a vector of nonnegative integers, a is a complex coefficient, and Fk{z) are non- 
trivial irreducible monic polynomials in z _1 . 

Defintion 2 (Conjugate Symmetry). A polynomial X(z) in z _1 is said to be conjugate 
symmetric if, for some vector k of positive integers and some 9 e [0, 2ir), 



X(z) = e^z^Xiz- 1 ). 

A conjugate symmetric polynomial may be reducible, irreducible, trivial, or nontrivial. If 
A(z) is an arbitrary polynomial in z _1 , then 



X(z) = A(z) ■ z" N A(z- 1 
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is conjugate symmetric. Any monomial az k is conjugate symmetric. 

The uniqueness of recovering a real- valued object from its Fourier magnitude or phase only 
data is discussed in [5] and can be easily generalized to the case of complex-valued objects. 

Proposition 4. Let /(n) G C(Af) be a finite array whose z-transform is irreducible up 
to a power of z _1 . If the Fourier transform G of g(n) G C(N) satisfies |G(e* 2pia; )| = 
|F(e i2 ™)|, G C, then 3 6 G [0, 2vr) and m such that either g = e ie f m+ or g = e ie J^l. 

Proposition 5. Let f G C{M) be a finite array whose z-transform has no nontrivial conju- 
gate symmetric factors. If g G C(N) satisfies lF(e 2ni ") = ^G(e 2m "),\/uj G C, then g = (3f 
for some real positive number (5. 

Proof. Consider the array h defined by 

h(n) = /(n) -k g(-n) 

whose z-transform is 

if (z) = F{z)G(z- 1 ). 
Since the phase of the Fourier transform of h(n) is equal to 

ZiJ(e 2 ™) = lF(e 2 ™") - ZG(e 2 ™), 

it follows that if lF(e 2 * iuj ) = lG(e 2niuj ), then lH(e 2 * ibJ ) = 0. Thus the Fourier transform 
of h is real-valued, implying that 

if (z) = H(z~ l ). 

Therefore, 

F(z)GWj = F^jG(z). (18) 
Multiplying both sides of f lTBj) by z~ N results in the following polynomial equation in z -1 : 

F(z)G(z- 1 )z- N = F(z 3T )G(z)z- N . (19) 

Since F(z) does not have trivial factors or nontrivial conjugate symmetric factors, we have 

F(z) = aJlF fc (z), (20) 

where ffc(z) are nontrivial irreducible non-conjugate symmetric monic polynomials in z _1 . 
Thus 

z- N F(z^) = a'z- m ' l[ F k (z), (21) 

k 
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where F^{z) are the nontrivial irreducible non-conjugate symmetric monic polynomials in 
z v of the form F^{z) = z _N+Pfe i 7 ) { .(z~ 1 ) for some vector of positive integers. 
Writing 

G(z) = 6z- n IjG < (z), (22) 

where Gg(z) are nontrivial irreducible monic polynomials in z _1 , we have 

z^G^ 1 ) = b'z- n ' H G t (z), (23) 

t 

where Gg(z) are the nontrivial irreducible monic polynomials in z _1 of the form Gg(z) = 
z~ t * +q£ Gg(z~ 1 ) for some vector of positive integers. 
Plugging flU, (J2TJ, (E2D and ((23]) in (TH yields 

a&'z""' J] F k (z) J] G,(z) = a'6z- m '-™ J] F fc (z) J] G,(z). (24) 

hi k e 

Each nontrivial irreducible factor F k (z) must be equal to some Fk'(z) or some G^(z). How- 
ever, if -Ffc(z) = -Ffc(z), then -Ffc(z) itself is conjugate symmetric. If, on the other hand, 
-Ffc(z) = -Ffc'(z) for some ^ k, i 7 fc(z)i 7 ) c /(z) = i 7 ^/ (z) Fy (z) becomes a conjugate symmetric 
factor. Both cases, however, are excluded by the assumption that the ^-transform of / does 
not have conjugate symmetric factors. Thus each F^{z) must be equal to G^(z) for some I' 
and Fiz) so that G(z) must be related by 

G(z) = Q(z) J] F k (z) = -Q(z)F(z). (25) 
k a 

However, G(z) and F(z) are both polynomials in z _1 , and since -F(z) contains no trivial 
factors, so Q{z) must be a polynomial of z -1 . Furthermore, plugging ( I2"5"j) in ffTSj) yields 

g( z ) = g(z^y. 

Therefore, Q(z) = /3 and the theorem follows by noting that /? must be positive if £F(oj) = 
ZG(w). 

□ 

We next show that the ^-transform of {A(n)/(n)} is almost surely irreducible up to a 
power z" 1 and not conjugate symmetric. 

Lemma 1. Let f G C(N) be a complex-valued array. Let {A(n)} be independent and con- 
tinuous random variables on S 1 . Then, V t ^ 0, the z-transform of (A/)t+ and (A/)t- is 
almost surely not conjugate symmetric. 
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Proof. Let 

; t+ (n) = A(t + n)/(t + n) 

whose ^-transform is 

F t+ (z) = A(t + n)/(t + n)z- n . (26) 

n 

F t+ (z) is conjugate symmetric if 

F t+ (z) = e ie z- k F t+ (z-i) (27) 
for some vector k of positive integers and some 9 £ [0, 2ir). Plugging (1261) in (|2"Tj) yields 
^A(t + n)/(t + n)z- n = e ie z~ u ^ A(t + n')/(t + n')z n ', 

n n' 

which implies 

A(t + n)/(t + n) = e ie A(t + k-n)/(t + k-n), Vn. (28) 

However, / is deterministic, and A(n) are independent and continuous random variables, 
so ( 128]) fails with probability one for any k. There are finitely many choices of k, so the 
z-transform of (A/)t+ is almost surely not conjugate symmetric. 

Similarly, the ^-transform of (A/) t _ is also almost surely not conjugate symmetric. □ 

Lemma 2. Let f £ C(J\f) be a complex-valued array of rank > 2. Let {A(n)} be independent 
and continuous random variables on S 1 . Then, the z-transform o/{A(n)/(n)} is irreducible 
up to a power of z _1 with probability one. 

For the proof of Lemma H] see Theorem 2 of [T] . 

Lemma 3. Let f and h be two complex-valued arrays. Let $ be the discrete Fourier operator 
such that = E k e " 2 ^' k /( k )- Then = implies that Z$/ = Z$/i ( _ t)+ . 

Proof. Note that 

$/t+M = e 27rit ' w $/(u;) 

which implies 

2?rt ■ u + Z$/(u>) (mod 2tt) = /.$h(u) 
by the assumption Z$/ t+ = Thus 

= Z.$h(w) - 27rt • u> (mod 2tt) 

which is equivalent to 

Z$/ = Z$/i(_ t)+ . 

□ 



18 



Let us now turn to the proof of Theorem HI 

Proof. Let / be the true image and h a fixed point of the ER iteration. Suppose that ti = Vfh 
satisfies the zero-padding condition. Then the following three equations hold: 

V h' = h (29) 
\$Ah'\ = |$A/| (30) 
l<$>Ati = Z<M/i (31) 

According to Lemma El the ^-transform of Af is irreducible up to a power of z -1 with 
probability one, so there exists some integer-valued vector m with — N < m < and some 
v G [0, 27r) such that 

ti = e i "A- 1 A m+ / in+ 

or 

In the case of h! = e^A _1 A m+ / m+ , the third equation in (I3"T|) becomes 

Ze^$A m+ / m+ = l^Ah. 

By Lemma [3], 

Ze^$A/ = Z$A ( _ m)+ / i( „ m)+ . (32) 

Lemma [H and El together with the assumption that /(0) 7^ imply that the z-transform 
of Af is an irreducible, nontrivial and non-conjugate symmetric polynomial of z" 1 with 
probability one. 

Next, we apply Proposition \5\ to (13"2"|) . Both Af and A(_ m ) + /i(_ m ) + are supported on a 
subset of {n I — N < n < N}. By Proposition we obtain 

ie w Af = A ( _ m)+ /i ( _ m _ )+ 

or equivalently 

h(n) = 7 e^ \ J f(n + m 
A(n) 

for some positive number /3. 

(a) If the true image /(•) is real- valued, then h = V h! is real- valued, which by the proof of 
Theorem [TJ (see Corrolary 1 of pQ) implies that v = 0, ir and m = or equivalently 

h = ±7/ (33) 

with probability one. Plugging (I3"3"|) in V Vfh = h yields 7 = 1 and thus h = ±/ with 
probability one. 
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(b) If / satisfies the sector condition of Theorem [21 then h = V a h' satisfies the same sector 
condition which by the proof of Theorem [2] (see Theorem 4 (i) of p]) implies that 
m = and 

h = je iv f (34) 

with probability at least 1 - \Af\(f3 - a)^ s ^(2n)^ s / 2 l Plugging flM} in V V f h = h 
yields 7 = 1 and thus h = e lv f. 

By the similar argument one reaches the same conclusion in the case of hi = 

□ 
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Fig. 3. (a) Recovery by 1651 ER iterations with UI and a = 4. (b) r(fk) versus 
k with r(f) « 5.49%. (c) Recovery by 1000 HIO + 103 ER with UI and a = 4. 
(d) r(ff C ) versus k with r(f) ~ 0.49%. (e) Recovery by 1587 ER steps with 
one low resolution RPI with a = 2. (f) r(fk) versus k with r(f) ~ 0.52% and 
e(j) ~ 2.51%. (g) Recovery by 33 HIO + 24 ER steps with low resolution 
RPI with (j — 2,. (h) r(f k ) versus k with r(f) w 0.05% and e{f) w 0.32%. 
(i) Recovery by 5512 ER steps with high resolution RPI with a = 1. (j) r(fk) 
versus k with r(f) w 0.19% and e(/) « 3.27%. (k) Recovery by 77 HIO + 
67 ER steps with high resolution RPI with a — 1. (1) r(//t) versus k with 
r(/) w 0.10% and e(/) w 1.39%. 
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(a) (b) (c) (d) 




(e) (f) (g) (h) 




(i) (j) (k) 



Fig. 4. (a) Recovery by 4140 ER iterations with UI with a = 4. (b) r(fk) versus 
k with r(f) « 14.71%. (c) Recovery by 1000 HIO + 421 ER steps with one UI 
with a = 4. (d) r(fk) versus k with r(f) ~ 3.94%. (e) Recovery by 460 ER with 
one low resolution RPI with a = 2. (f) r(fk) versus k with r(f) ~ 0.03% and 
e(/) ~ 0.09%. (g) Recovery by 103 HIO + 11 ER steps with one low resolution 
RPI with a = 2. (h) r(f k ) versus k with r(f) » 0.03% and e(/) w 0.12%. (i) 
Recovery by 966 ER steps with one high resolution RPI with a = 1. (j) 
versus k with r(f) w 0.06% and e(/) w 0.40%. (k) Recovery by 94 HIO + 
16 ER steps with one high resolution RPI with a — 1. (1) r(/fc) versus k with 
r(f) w 0.04% and e(/) w 0.19%. 
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relative error versus oversampling rate 



relative error versus oversampling rate 




high resolution & noise free 
low resolution & noise free 
high resolution & 5% gaussian nosie 
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(c) 



Fig. 5. (a) Relative error with one RPI for nonnegative- valued Phantom; (b) 
Relative error with one RPI for complex-valued Phantom with phases ran- 
domly distributed in [0, 7r/2]; (c) Relative error by 200 HIO +300 ER with one 
RPI and UI for complex-valued Phantom with phases randomly distributed in 
[0, 2tt} with one UI and one RPI of high resolution (block size: 1 x 1) or low 
resolution (block size: 4x4). 
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(a) (b) (c) (d) 




(e) (f) (g) (h) 




(i) (j) (k) (1) 



Fig. 6. Phasing with a — 2 and one high resolution RPI: (a) Recovery by 18 
HIO +10 ER with 5% Gaussian noise; (b) r(f k ) versus k with r(f) « 2.62% 
and e(/) « 4.20%; (c) Recovery by 19 HIO +10 ER with 5% Gaussian noise, 
(d) r(/jfc) versus k with r(f) ~ 2.85% and e(f) ~ 3.51%; (e) Recovery by 16 
HIO +10 ER with 5% Poisson noise; (f) r(/fc) versus k with r(f) « 3.71% and 
e(/) « 5.89%; (g) Recovery by 17 HIO +10 ER with 5% Poisson noise; (h) 
r(f k ) versus k with r(f) m 4.05% and e(f) « 4.84%; (i) Recovery by 14 HIO 
+ 10 ER with 5% illuminator noise; (j) r(fk) versus k with r(f) ~ 5.28% and 
e(f) w 7.75%; (k) Recovery by 16 HIO +10 ER with 5% illuminator noise; (1) 
r(fk) versus k with r(f) w 5.48% and e(/) w 6.35%. 



25 




(a) (b) (c) (d) 




(e) (f) (g) (h) 




(i) (j) (k) (1) 



Fig. 7. Phasing with a = 2 and one low (40 x 40) resolution RPI: (a) Recovery 
by 27 HIO +10 ER with 5% Gaussian noise; (b) r(fk) versus k with r(f) ~ 
2.50% and e(f) « 7.37%; (c) Recovery by 35 HIO +10 ER with 5% Gaussian 
noise, (d) r(fk) versus k with r(f) ~ 2.85% and e(/) ~ 4.18%; (e) Recovery by 
22 HIO +10 ER with 5% Poisson noise; (f) r(f k ) versus k with r(/) « 3.77% 
and e(f) « 6.27%; (g) Recovery by 100 HIO +10 ER with 5% Poisson noise; 
(h) r(fk) versus k with r(f) « 4.24% and e(f) « 5.09%; (i) Recovery by 20 
HIO +10 ER with 5% illuminator noise; (j) r(fk) versus k with r(f) ~ 4.00% 
and e(f) w 13.14%; (k) Recovery by 34 HIO +10 ER with 5% illuminator 
noise; (1) r(f k ) versus k with r(f) ~ 5.48% and e(/) ~ 9.46%. 
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relative error versus noise level 
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Fig. 8. (a) Relative error for nonnegative-valued Phantom and a = 2 (b) 
Relative error for complex-valued Phantom with phases randomly distributed 
in [0,7r/2] and a = 4; (c) Relative error for complex- valued Phantom with 
phases randomly distributed in [0, 2tt] and a = 3. 
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